Bienvenid@s a la tercera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en el diseño de experimentos, test de hipótesis y regresión lineal. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Enlaces a videos de las clases:
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 3 son las siguientes:
# Manipulación de estructuras
library(tidyverse)
library(dplyr)
library(tidyr)
# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
# Manipulación de varios plots en una imagen.
library(gridExtra)
En clases se han visto diferentes tipos de test de hipótesis para demostrar una proposición sobre algún parámetro. Uno de los test vistos en clases es el Z-Test, el cual su distribución del test estadístico bajo la hipótesis nula se puede aproximar a una Gaussina. Para la aplicación de este test, resaltan los siguientes puntos:
Para calcular la significancia estadística al igual que con otros métodos esta se debe calcular como:
Luego, dependiendo del objetivo del test tenemos las metodologías one-sample y two-sample. Utilizaremos One-Sample cuando nuestro objetivo es comparar la media de una muestra con la media de la población. El Z-score del One-Sample se define como:
\[Z-score_{One-Sample} = \dfrac{\bar x - \mu}{\dfrac{\sigma}{\sqrt n}}\] Donde \(\bar x\) es la media de la muestra, \(\mu\) es la media de la población, \(\sigma\) es la desviación estándar de la población y \(n\) es el tamaño de la muestra.
Por otro lado, se utiliza Two-Sample cuando queremos comparar la media de dos muestras. El Z-score de Two-Sample se define con la ecuación:
\[Z-score_{Two-Sample} = \dfrac{(\bar x_1 - \bar x_2) - (\mu_1 - \mu_2)}{\sqrt{(\dfrac{\sigma_1}{\sqrt n_1}+\dfrac{\sigma_2}{ n_2}})}\] Donde \((\bar x_2 - \bar x_1)\) es la diferencia de las medias de la muestra, \((\mu_1 - \mu_2)\) la diferencia de las medias de la población, \(\sigma_{1,2}\) la desviación estándar de la población y \(n_{1,2}\) el tamaño de las muestras.
En la práctica aparece la necesidad de testear múltiples hipótesis (por ejemplo en biología se pueden utilizar múltiples grupos de control o querer estudiar múltiples resultados de un mismo experimento), de esta forma la primera idea es testear individualmente cada una de las hipótesis, el problema de este enfoque es que la probabilidad de que se obtenga al menos un resultado significante crece rápidamente (con un nivel de significancia \(\alpha = 0.05\) y \(20\) test ya se alcanza una probabilidad de \(64\%\) de tener resultados significantes por azar).
Una forma de corregir los inconvenientes del método anterior es utilizar el método de Bonferroni correction quien propone cambiar \(\alpha\) por \(\alpha/m\) (donde \(m\) es la cantidad de test de hipotesis realizados), esto resulta que las probabilidades de rechazar por error se mantengan bajas. De esta forma los p-valores obtenidos en un test de hipótesis y al utilizar Bonferroni correction, quedan dados por el producto de un \(p-valor_{i}\) y la cantidad de test realizados: \(\text{p-valor}_{i}*m\).
El objetivo de esta pregunta es programar la potencia de un test de hipótesis y observar como se comportan las la hipótesis nula v/s la alternativa para un Z-test. Con el desarrollo de este ejercicio, podrán visualizar las diferentes partes que conforman a un test de hipótesis, identificar que es el p-valor y evidenciar como varia la potencia de un test one-sample y two-sample al variar \(\alpha\).
Para recordar; sabemos que en estadística el concepto de potencia viene dado por:
\[Power = 1 - \beta\]
Donde \(\beta\) es la probabilidad de obtener un error de tipo II. Con esto, la potencia estadística viene a representar la probabilidad de rechazar la hipótesis nula cuando esta es falsa. O sea, la potencia de una prueba es la probabilidad de encontrar un resultado positivo dado que este existe. Una de las formas de representar la potencia de un test es a través del siguiente gráfico:
Del gráfico, es posible visualizar que a medida que aumenta la diferencia en la media de la población, se obtienen mayores valores de potencia estadística.
Recordada que es la potencia de un test de hipótesis, a continuación, usted deberá programar una función que sea capaz de obtener la potencia de un Z-test one-sample y two-sample. Para esto por favor considere los siguientes puntos:
function(n1=NULL, sigma1=0.5,
n2=NULL,sigma2=0.5, mu.Ha=0 ,
mu.True=0, alfa=0.05)
De los argumentos, tendremos que: \(n1\) representa la cantidad de datos para la muestra 1, \(sigma1\) es la desviación estándar de la muestra 1, \(n2\) la cantidad de datos para la muestra 2, \(sigma2\) la desviación estándar para la muestra 2, \(mu.Ha\) el mu del test de hipótesis y \(mu.True\) la media de la población real. Notar que la presencia de una segunda muestra solo es para el caso two-sample, para el caso one-sample el argumento de entrada \(n2\) debería ser nulo.
Codificada la función realice los siguientes experimentos:
\[ n1=16, sigma1=16, mu.Ha=100 , mu.True=Variar, alfa=0.05 \] \[ n1=16, sigma1=16, mu.Ha=100 , mu.True= Variar, alfa=0.01 \] \[ n1=16, sigma1=16, mu.Ha=100 , mu.True= Variar, alfa=0.1 \]
Se le recomienda que la variación se realice a través de un
for y grafique las curvas dentro de un mismo gráfico para
observar potenciales diferencias entre ellas.
Diseñe un experimento one-sample y visualice cómo se comportan las distribuciones normales de la hipótesis nula y la hipótesis alternativa al variar \(\alpha\).
Diseñe un experimento Two-sample y visualice cómo se comportan las distribuciones normales de la hipótesis nula y la hipótesis alternativa al variar \(\alpha\).
Para el diseño de experimentos y/o comprobación de sus métodos puede serles útiles (no hay problema si decide utilizar los mismos ejemplos):
Respuesta
# Power Function, El esqueleto posee como ejemplo como obtener la potencia de un z-test one-sample.
# Si utiliza este esqueleto deberá comentar la función que cumple cada una de las partes entregadas
power.z.test <- function(n1=NULL, sigma1=0.5,
n2=NULL,sigma2=0.5, mu.Ha=0 ,
mu.True=0, alfa=0.05){
if(){
#
Z = qnorm(1-alfa)
denominador = sigma1/sqrt(n1)
X_bar = Z*denominador + mu.Ha
numerador = X_bar - mu.True
Z = numerador/denominador
Power = 1 - pnorm(Z)
#
min_lim = min(rnorm(1000, mean=mu.Ha, sd=denominador)) -
round(min(rnorm(1000, mean=mu.Ha, sd=denominador)))%%10
max_lim = max(rnorm(1000, mean=mu.True, sd=denominador)) +
round(max(rnorm(1000, mean=mu.True, sd=denominador)))%%10
#
plot <- ggplot(data.frame(x = c(min_lim, max_lim)), aes(x)) +
stat_function(fun = dnorm, args = list(mean = mu.Ha, sd = denominador),
col='red') +
stat_function(fun = dnorm, args = list(mean = mu.True, sd = denominador),
col='blue') +
stat_function(fun = dnorm, args = list(mean = mu.True, sd = denominador),
xlim = c(X_bar,max_lim), geom = "area", fill='red') +
geom_vline(xintercept = X_bar, linetype="dotted", size=1) +
annotate(x=X_bar, y=+Inf,label="alpha", vjust=2, geom="label") +
theme_minimal() +
ggtitle("H0 vs Ha") +
xlab(expression(bar(X))) + ylab("Density")
}
}
if(){
#
#
if(plot.hypothesis){
#
#
plot = ...
}
}
# Como R no permite retornar dos salidas usamos una lista
# Los resultados se llaman con $plot o $power
return(list(plot=plot,power=Power))
}
# Plot de gráfico de potencia
# Experimentos
Esta pregunta tiene como objetivo comprender como funciona un test de hipótesis y como deberíamos abordar la realización de múltiples test de hipótesis con datos reales.
La pregunta deberá ser desarrollada utilizando el dataset
marketing_campaign.csv. Con esto, deberá programar un
Z-test, con el cual estudiará a través de experimentos el
Income de personas con los grados académicos
Graduation, Master y PhD. Para
realizar esto considere la elaboración de los siguientes puntos de forma
secuencial:
Graduation, Master y PhD por
separado.
Por ejemplo para el caso de Graduation pueden generar estructuras de la siguiente forma:
| ID | Graduation |
|---|---|
| 5524 | 58138 |
| 2174 | 46344 |
| 4141 | 71613 |
| 6182 | 26646 |
| 965 | 55635 |
| … | … |
Donde los valores en la fila de Graduation representan los sueldos de las diferentes personas que conforman el dataset. Un punto importante a considerar es que los datos para los diferentes grados académicos poseen diferentes numero de datos (no se asusten por esto).
Programar el método Z-test con la metodología one sample y two sample, obteniendo los p-valores a través de las alternativas one-sided y two-sided. Para el caso de one-sided, cree una función capaz de obtener la cola menor y mayor de la gaussiana.
El calculo de las diferentes alternativas para calcular los p-valores deberá ser un argumento de su función, donde señalando ‘menor’,‘mayor’ (para los casos one-sided) y ‘two-sided’ deberá obtener el valor pertinente para cada caso.
Genere una función que permita realizar solo múltiples test del
tipo two-sample y aplique bonferroni correction a los p-valores
obtenidos. Notar que los múltiples test deberá realizar la comparación
entre todos los elementos de entrada, por ejemplo si deseamos comparar
los ingresos de Graduation, Master y
PhD, se deberían comparar los ingresos de
Graduation v/s Master, Graduation
v/s PhD y Master y PhD
Codificada las funciones, realice los siguientes experimentos con su función de test de hipótesis:
Compruebe si la media de los ingresos para la variable
Graduation es similar a 52000. Señale formalmente este
experimento y obtenga los p-valores para las alternativas one-sided y
two-sided.
Compruebe si la diferencia entre los ingresos de las personas con
el grado académico Graduation es cercana a cero en relación
a la recibida por los Master y PhD. Para este
punto utilice la función que le permite realizar múltiples test del tipo
two-sample.
Para los diferentes experimentos considere que la desviación estandar
de la población para los diferentes income son los
siguientes:
\[\sigma_{Graduation} = 28180\] \[\sigma_{Master} = 20160\] \[\sigma_{PhD} = 20615\]
Respuesta:
Importación y primera visualización de los datos
df = read.csv('marketing_campaign.csv', sep='\t')
graduation <- df[df$Education == "Graduation", "Income"]
master <- df[df$Education == "Master", "Income"]
phd <- df[df$Education == "PhD", "Income"]
head(df)
Definición de la función z_test
# Implementación de Z-test one-sided y two-sided
z_test <- function(data1 = NULL, data1_name = NULL, sigma1 = 0.5,
data2 = NULL, data2_name = NULL, sigma2 = 0.5,
mu1 = 0, mu2 = 0, verbose = TRUE,
test_type = c("one-sided", "two-sided")) {
# Condiciones para evitar errores en la función
length_test_type_condition <- length(test_type) == 1
name_test_type_condition <- (test_type %in% c("menor", "mayor", "two-sided"))
if (!(length_test_type_condition && name_test_type_condition)) {
print("Por favor escoge un tipo de Test: ´mayor´, ´menor´ o ´two-sided´ ")
return()
} else if (is.null(data1)) {
print("Por favor ingresa algún valor en data1")
return()
}
# Nombre data1
if (is.null(data1_name))
data_deparse <- deparse(substitute(data1))
else
data_deparse <- data1_name
# Z-Score
if (is.null(data2)) {
mean_data <- mean(data1, na.rm = TRUE)
z_score <- (mean_data - mu1) / (sigma1 * sqrt(length(data1)))
output <- "One"
} else {
mean_data1 <- mean(data1, na.rm = TRUE)
mean_data2 <- mean(data2, na.rm = TRUE)
n_1 <- length(data1)
n_2 <- length(data2)
mean_diff <- (mean_data1 - mean_data2) - (mu1 - mu2)
z_score <- mean_diff / sqrt(sigma1^2 / n_1 + sigma2^2 / n_2)
# Nombre data2
output <- "Two"
if (is.null(data2_name))
data_deparse <- paste(data_deparse, "y", deparse(substitute(data2)))
else
data_deparse <- paste(data_deparse, "y", data2_name)
}
# P-Value
if (test_type == "menor")
p_value <- pnorm(z_score)
else if (test_type == "mayor")
p_value <- 1 - pnorm(z_score)
else if (test_type == "two-sided")
p_value <- 2 * pnorm(-abs(z_score))
# Texto de Salida
if (verbose) {
cat("\n\t", output, "-sample Z-Test:\n\nData analizada: ",
data_deparse, "\nZ = ", z_score,
" P-value = ", p_value,
"\n\n", sep = ""
)
}
return(p_value)
}
Implementación de la función z_test_multiple_testing, la
cual llama a la función anterior z_test
z_test_multiple_testing <- function(data = NULL, sigma = NULL,
test_type = c("one-sided", "two-sided"),
verbose = TRUE) {
# Condiciones para evitar errores en la función
length_test_type_condition <- length(test_type) == 1
name_test_type_condition <- (test_type %in% c("menor", "mayor", "two-sided"))
nullity_condition <- is.null(data) || is.null(sigma)
dimension_condition <- !nullity_condition && (length(data) == length(sigma))
if (!(length_test_type_condition && name_test_type_condition)) {
print("Por favor escoge un tipo de Test: ´mayor´, ´menor´ o ´two-sided´ ")
return()
} else if (!dimension_condition) {
cat(
"Cantidad incompatible de datos y sigmas\n",
"length(data) = ", length(data), "\n",
"length(sigma) = ", length(sigma), "\n\n", sep = ""
)
}
# Nombres de los datos estudiados
data_names <- sapply(substitute(data), deparse)[-1]
m <- 0
# Cálculo del p-value para cada par de datos
indexes <- seq_along(data)
p_value_list <- list()
for (i in indexes) {
for (j in indexes) {
if (j <= i)
next
var1 <- data_names[i]
var2 <- data_names[j]
p_value <- z_test(
data1 = data[[i]], data2 = data[[j]],
data1_name = var1, data2_name = var2,
sigma1 = sigma[i], sigma2 = sigma[j],
test_type = test_type, verbose = verbose
)
key <- paste(var1, "_", var2, sep = "")
p_value_list[key] <- p_value
m <- m + 1
}
}
# Corrección de Bonferroni
p_value_list <- lapply(p_value_list, FUN = function(x) x * m)
return(p_value_list)
}
Para realizar el primer experimento, el cual consiste en determinar
si la media de la variable Graduation es similar a 52000,
realizaremos 3 tests de hipótesis, siempre considerando que \(\mu_{Grad}\) es la media de la variable
Graduation, \(\mu_0 =
52000\) y \(\alpha = 0.05\).
Para el primer test, se usarán las siguientes hipótesis: \[H_0 : \mu_{Grad} = 52000\] \[H_a : \mu_{Grad} \neq 52000\]
Para este test, ya que la hipótesis nula \(H_0\) tiene una igualdad, se utilizará un test de dos colas.
p_value <- z_test(data1 = graduation, mu1 = 52000,
test_type = "two-sided", sigma1 <- 28180,
data1_name = "graduation")
##
## One-sample Z-Test:
##
## Data analizada: graduation
## Z = 0.0007614736 P-value = 0.9993924
output <- if (p_value <= 0.05) "Se rechaza" else "No se rechaza"
cat(output, "la hipótesis nula a un nivel de 5%")
## No se rechaza la hipótesis nula a un nivel de 5%
Para el segundo test, se usarán las siguientes hipótesis: \[H_0 : \mu_{Grad} \geq 52000\] \[H_a : \mu_{Grad} < 52000\]
Para este test, ya que la hipótesis nula \(H_0\) tiene un \(\geq\), se utilizará un test de una cola izquierda.
p_value <- z_test(data1 = graduation, mu1 = 52000,
test_type = "menor", sigma1 <- 28180,
data1_name = "graduation")
##
## One-sample Z-Test:
##
## Data analizada: graduation
## Z = 0.0007614736 P-value = 0.5003038
output <- if (p_value <= 0.05) "Se rechaza" else "No se rechaza"
cat(output, "la hipótesis nula a un nivel de 5%")
## No se rechaza la hipótesis nula a un nivel de 5%
Finalmente, para el tercer test, se usarán las siguientes hipótesis: \[H_0 : \mu_{Grad} \leq 52000\] \[H_a : \mu_{Grad} > 52000\]
Para este test, ya que la hipótesis nula \(H_0\) tiene un \(\leq\), se utilizará un test de una cola derecha.
p_value <- z_test(data1 = graduation, mu1 = 52000,
test_type = "mayor", sigma1 <- 28180,
data1_name = "graduation")
##
## One-sample Z-Test:
##
## Data analizada: graduation
## Z = 0.0007614736 P-value = 0.4996962
output <- if (p_value <= 0.05) "Se rechaza" else "No se rechaza"
cat(output, "la hipótesis nula a un nivel de 5%")
## No se rechaza la hipótesis nula a un nivel de 5%
De estos 3 tests se puede concluir que, como la hipótesis nula nunca
se rechazó, entonces se acepta como la correcta. Es decir, la media de
la variable Graduation es similar a 52000.
Para estudiar la diferencia entre los ingresos de las personas dependiendo del grado académico, se va a realizar un test múltiple sobre los siguientes conjuntos de hipótesis:
Sean \(\mu_{grad}\), \(\mu_{master}\) y \(\mu_{phd}\) las medias de los datos y \(\alpha = 0.05\)
\[\begin{cases} H_1: \mu_{grad} = \mu_{master} \\ \bar{H_1}: \mu_{grad} \neq \mu_{master} \end{cases}\]
\[\begin{cases} H_2: \mu_{grad} = \mu_{phd} \\ \bar{H_2}: \mu_{grad} \neq \mu_{phd} \end{cases}\]
\[\begin{cases} H_3: \mu_{master} = \mu_{phd} \\ \bar{H_3}: \mu_{master} \neq \mu_{phd} \end{cases}\]
p_values_list <- z_test_multiple_testing(
data = list(graduation, master, phd),
sigma = c(28180, 20160, 20615),
test_type = "two-sided"
)
##
## Two-sample Z-Test:
##
## Data analizada: graduation y master
## Z = -0.1468296 P-value = 0.8832666
##
##
## Two-sample Z-Test:
##
## Data analizada: graduation y phd
## Z = -2.725542 P-value = 0.0064196
##
##
## Two-sample Z-Test:
##
## Data analizada: master y phd
## Z = -2.298014 P-value = 0.021561
grad_master_text <- if (p_values_list["graduation_master"] <= 0.05) "Se" else "No se"
grad_phd_text <- if (p_values_list["graduation_phd"] <= 0.05) "Se" else "No se"
master_phd_text <- if (p_values_list["master_phd"] <= 0.05) "Se" else "No se"
cat("",
grad_master_text, "rechaza H_1\n",
grad_phd_text, "rechaza H_2\n",
master_phd_text, "rechaza H_3\n"
)
## No se rechaza H_1
## Se rechaza H_2
## No se rechaza H_3
Podemos notar que, dadas las hipótesis rechazadas y las que no, la
media de los ingresos de las personas con grado Graduation
es muy cercana a la media de las personas con grado Master,
pero lejana a la media de las personas con grado PhD. No
obstante, es importante recalcar que la media de los grado
Master es cercana también a la media de los grado
PhD, lo cual nos indica que la media de ingresos de
Master está entre ambas, y por ello a una “distancia”
razonable a las otras dos.
El objetivo de este problema es estudiar como realizar múltiples test
de hipótesis simultáneamente. Para esto en primer lugar se estudiara el
método “intuitivo”, donde veremos sus limitantes y se comparará con el
método llamado Bonferroni correction, posteriormente se
realizará un estudio practico con el dataset
ratones.csv.
Un investigador se ha colocado en contacto con ustedes señalándoles que realiza diariamente test de hipótesis entre las muestras que toma día a día en su laboratorio. Con esto, al investigador le urge saber si realizar multiples test de hipótesis sin una corrección podría afectar la toma de decisiones. Para comprobar esto, les solicita comprobar matemáticamente como se comporta la probabilidad de obtener al menos un resultado significativos al azar de sus experimentos diarios. Para esto, les señala que la la probabilidad de obtener un experimento por azar puede ser simulado a través de los casos exitosos de una binomial (valores mayores a cero), donde el numero de observaciones son la cantidad de experimentos (\(m\)) y la probabilidad queda dada por \(\alpha\) del test.
A continuación, se entregan unas indicaciones mas especificas para desarrollar la pregunta:
data <- read.csv("ratones.csv",sep= ";", stringsAsFactors = T)
head(data)
Respuesta Aquí:
probEmpirica <- function(alpha,m){
n <- 10000 # Cantidad de veces que se va a repetir el experimento para estimar la probabilidad, pueden cambiar este valor si lo desean
res <-... #Resultados de los experimentos
# Puede agergar todo el codigo que estime conveniente para calcular la probabilidad empirica
...
...
#
prob <-... # Probabilidad empirica
return(prob)
}
El objetivo de la siguiente pregunta es aplicar los conceptos de regresión lineal vistos en clases para implementar desde cero un función capaz de realizar una regresión simple y múltiple.
Para este problema, ustedes deberán estudiar el comportamiento de los
clientes de un holding de salud. Para esto, se les hace entrega del
dataset insurance.csv para que estudien la creación de un
modelo lineal con sus datos. Antes de comenzar a trabajar, se señalan
las diferentes variables que componen al dataset:
Es importante que considere que cada una de las filas representa un cliente distinto para el holding.
Dentro del estudio, el holding de salud le solicita estudiar los
comportamientos de los clientes fumadores y no fumadores, por lo que se
le aconseja separar el dataframe original en fumadores y no fumadores.
En el estudio, realicen un modelo lineal que tiene como variable de
respuesta a charges y los datos que mejor se correlacionan
para los clientes fumadores y no fumadores. Para esto, deberán realizar
las siguientes actividades:
charges.lm(), ya que esto le servirá para observar si la
regresión lineal es significativa.Nota: No esta permitido utilizar comandos que obtengan los valores solicitados directamente a menos que se le permita en la pregunta.
Desarrollo:
El modelo lineal está implementado usando la formulación matemática de los coefficientes: \[\hat{\beta_1} = \frac{\sum\limits_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum\limits_{i=1}^n (x_i - \bar{x})^2}\\ \hat{\beta_0} = \bar{y} - \hat{\beta_1}\bar{x}\] donde \(\hat{\beta_0}\) es el intercepto y \(\hat{\beta_1}\) es la pendiente de la recta.
linear_regression <- function(data, target) {
x_mean <- mean(data)
y_mean <- mean(target)
numerator <- sum((data - x_mean) * (target - y_mean))
denominator <- sum((data - x_mean)^2)
beta_1 <- numerator / denominator
beta_0 <- y_mean - beta_1 * x_mean
list("intercept" = beta_0, "slope" = beta_1)
}
Luego, se pasan las variables categóricas sex y
region a variables dummy, para después separar los datos
entre fumadores y no fumadores
data <- read.csv("insurance.csv")
# One Hot Encoding
data$female <- ifelse(data$sex == "female", 1, 0)
data$male <- ifelse(data$sex == "male", 1, 0)
data$southwest <- ifelse(data$region == "southwest", 1, 0)
data$southeast <- ifelse(data$region == "southeast", 1, 0)
data$northwest <- ifelse(data$region == "northwest", 1, 0)
data$northeast <- ifelse(data$region == "northeast", 1, 0)
data <- data[, !names(data) %in% c("sex", "region")]
# Separación de datos en 'smoker' y 'non_smoker'
smoker <- data[data$smoker == "yes", !names(data) == "smoker"]
non_smoker <- data[data$smoker == "no", !names(data) == "smoker"]
Una vez separados los datos, se obtiene la variable con mayor
correlación a la variable charges, y se preparan los datos
para ser ingresados al modelo lineal
get_max_correlation_column <- function(df, column) {
cor_df <- as.data.frame(as.table(cor(df)))
cor_df <- cor_df[cor_df$Var1 == column, ]
cols <- c("Var1", "Var2")
cor_df[cols] <- lapply(cor_df[cols], as.character)
max_cor <- max(abs(cor_df[cor_df$Freq != 1, ]$Freq))
variable <- cor_df[abs(cor_df$Freq) == max_cor, ]$Var2
variable
}
smoker_variable <- get_max_correlation_column(smoker, "charges")
x_smoker <- smoker[, smoker_variable]
y_smoker <- smoker$charges
non_smoker_variable <- get_max_correlation_column(non_smoker, "charges")
x_non_smoker <- non_smoker[, non_smoker_variable]
y_non_smoker <- non_smoker$charges
Con los datos listos, los ingresamos al modelo lineal y graficamos los resultados
# Regresión Lineal para los datos 'smoker'
reg_smoker <- linear_regression(x_smoker, y_smoker)
slope_smoker <- reg_smoker$slope
intercept_smoker <- reg_smoker$intercept
plot(x_smoker, y_smoker, ylab = "charges", xlab = smoker_variable)
abline(a = intercept_smoker, b = slope_smoker, col = "red")
title("Linear Regression for 'smokers'")
# Regresión Lineal para los datos 'non_smoker'
reg_non_smoker <- linear_regression(x_non_smoker, y_non_smoker)
slope_non_smoker <- reg_non_smoker$slope
intercept_non_smoker <- reg_non_smoker$intercept
plot(x_non_smoker, y_non_smoker, ylab = "charges", xlab = non_smoker_variable)
abline(a = intercept_non_smoker, b = slope_non_smoker, col = "red")
title("Linear Regression for 'non_smokers'")
Veamos el coeficiente \(R^2\) para cada caso
# Coeficiente R2 para 'smokers'
y_hat_smoker <- intercept_smoker + slope_smoker * x_smoker
SST_smoker <- sum((y_smoker - mean(y_smoker))^2)
SSE_smoker <- sum((y_smoker - y_hat_smoker)^2)
R2_smoker <- 1 - (SSE_smoker / SST_smoker)
print(sprintf("R2 for 'smoker': %.3f", R2_smoker))
## [1] "R2 for 'smoker': 0.650"
# Coeficiente R2 para 'non_smokers'
y_hat_non_smoker <- intercept_non_smoker + slope_non_smoker * x_non_smoker
SST_non_smoker <- sum((y_non_smoker - mean(y_non_smoker))^2)
SSE_non_smoker <- sum((y_non_smoker - y_hat_non_smoker)^2)
R2_non_smoker <- 1 - (SSE_non_smoker / SST_non_smoker)
print(sprintf("R2 for 'non_smoker': %.3f", R2_non_smoker))
## [1] "R2 for 'non_smoker': 0.394"
Ahora el coeficiente \(\bar{R^2}\), o \(R^2\) ajustado
# Coeficiente R2 ajustado para 'smokers'
n_smoker <- length(y_smoker)
R2_adj_smoker <- 1 - (1 - R2_smoker) * (n_smoker - 1) / (n_smoker - 2)
print(sprintf("R2 ajustado for 'smoker': %.3f", R2_adj_smoker))
## [1] "R2 ajustado for 'smoker': 0.649"
# Coeficiente R2 ajustado para 'non_smokers'
n_non_smoker <- length(y_non_smoker)
R2_adj_non_smoker <- 1 - (1 - R2_non_smoker) * (n_non_smoker - 1) / (n_non_smoker - 2)
print(sprintf("R2 ajustado for 'non_smoker': %.3f", R2_adj_non_smoker))
## [1] "R2 ajustado for 'non_smoker': 0.394"
Ahora, se implementará el modelo lineal con dos variables. Para esto,
se modificó la función que encuentra las variables con mayor correlación
con la variable charges y la función que encuentra el
modelo lineal, ya que ahora para encontrar el vector \(\hat{\beta}\) se usó la siguiente ecuación:
\[\hat{\beta} = (X^TX)^{-1}XY\]
get_max_correlation_column <- function(df, column) {
cor_df <- as.data.frame(as.table(cor(df)))
cor_df <- cor_df[cor_df$Var1 == column, ]
cols <- c("Var1", "Var2")
cor_df[cols] <- lapply(cor_df[cols], as.character)
max_cor_1 <- max(abs(cor_df[cor_df$Freq != 1, ]$Freq))
next_max <- abs(cor_df$Freq) < abs(max_cor_1)
max_cor_2 <- max(abs(cor_df[next_max, ]$Freq))
variable_1 <- cor_df[abs(cor_df$Freq) == max_cor_1, ]$Var2
variable_2 <- cor_df[abs(cor_df$Freq) == max_cor_2, ]$Var2
list("var1" = variable_1, "var2" = variable_2)
}
linear_regression_multivar <- function(data, target) {
matrix_data <- data.matrix(data)
colnames(matrix_data) <- NULL
beta <- solve(t(matrix_data) %*% matrix_data) %*% t(matrix_data) %*% target
list("beta_0" = beta[1], "beta_1" = beta[2], "beta_2" = beta[3])
}
Luego se entrenó el modelo
smoker_variables <- get_max_correlation_column(smoker, "charges")
non_smoker_variables <- get_max_correlation_column(non_smoker, "charges")
# Entrenamiento para 'smoker'
x_smoker <- data.matrix(smoker[, names(smoker) %in% smoker_variables])
# Se agrega una columna de 1's al principio
x_smoker <- cbind(rep(1, nrow(x_smoker)), x_smoker)
y_smoker <- smoker[, "charges"]
reg_smoker <- linear_regression_multivar(x_smoker, y_smoker)
beta_smoker <- unlist(reg_smoker)
y_hat_smoker <- x_smoker %*% beta_smoker
SST_smoker <- sum((y_smoker - mean(y_smoker))^2)
SSE_smoker <- sum((y_smoker - y_hat_smoker)^2)
R2_smoker <- 1 - (SSE_smoker / SST_smoker)
n_smoker <- length(y_smoker)
R2_adj_smoker <- 1 - (1 - R2_smoker) * (n_smoker - 1) / (n_smoker - 3)
writeLines(sprintf("R2 for 'smoker': %.3f\nR2 ajustado for 'smoker': %.3f", R2_smoker, R2_adj_smoker))
## R2 for 'smoker': 0.753
## R2 ajustado for 'smoker': 0.751
# Entrenamiento para 'non_smoker'
x_non_smoker <- data.matrix(non_smoker[, names(non_smoker) %in% non_smoker_variables])
# Se agrega una columna de 1's al principio
x_non_smoker <- cbind(rep(1, nrow(x_non_smoker)), x_non_smoker)
y_non_smoker <- non_smoker[, "charges"]
reg_non_smoker <- linear_regression_multivar(x_non_smoker, y_non_smoker)
beta_non_smoker <- unlist(reg_non_smoker)
y_hat_non_smoker <- x_non_smoker %*% beta_non_smoker
n_non_smoker <- length(y_non_smoker)
R2_adj_non_smoker <- 1 - (1 - R2_non_smoker) * (n_non_smoker - 1) / (n_non_smoker - 3)
writeLines(sprintf("R2 for 'non_smoker': %.3f\nR2 ajustado for 'non_smoker': %.3f", R2_non_smoker, R2_adj_non_smoker))
## R2 for 'non_smoker': 0.394
## R2 ajustado for 'non_smoker': 0.393
A work by CC6104